Denoising method and system for preserving clinically significant structures in reconstructed images using adaptively weighted anisotropic diffusion filter

ABSTRACT

Embodiments and processes of computer tomography perform tasks associated with denoising a reconstructed image using an anisotropic diffusion filter and adaptively weighting an iterative instance of the diffused image based upon the product of a weight value and a difference between the iterative instance of the diffused image and the original image. In general, the adaptive weighting is a negative feedback in the iterative steps.

CROSS REFERENCE TO RELATED APPLICATION

This application is a Divisional Application of U.S. application Ser.No. 12/913,796, filed on Oct. 28, 2010, the entire contents of which areincorporated herein by reference.

FIELD

Embodiments described herein generally relate to denoising in imagingsystems and methods.

BACKGROUND

As computed tomography (CT) achieves higher resolution, noise has been asignificant problem from the inception of CT imaging. A major concern indeveloping a denoising technique is to preserve clinically significantstructures of interest. In this regard, although low-pass filteringtechniques are well-known effective denoising filters, the processedimages tend to lose clinically significant structures. That is, low-passfiltering techniques generally remove high frequencies from the image inorder to reduce noise for improving the detectability of large objectsin the image. On the other hand, the low-pass filters also reduce theintrinsic resolution of the image as they smooth edges and consequentlydecrease the detectability of small structures.

Some improvements have been made to the low-pass filtering techniques.For example, some prior art low-pass filtering methods average voxels inlocal neighborhoods that are within an intensity range depending on theestimated standard deviation of noise throughout the image. Other priorart low-pass filtering techniques are based upon the iterativelyestimated standard deviation of noise in a particular local position ofCT volume for performing a weighted average between neighboring voxels.These noise reduction methods are not particularly aimed at preservingspecific structures of interest.

In contrast, a recent CT denoising approach involves anisotropicdiffusion as a denoising filter whose diffusion coefficient is madeimage dependent. Certain image structures can be directionally smoothedwith anisotropic diffusion. Since exemplary CT noise reduction methodsusing anisotropic diffusion improve only certain predeterminedstructures such as elongated-shaped vessels, these methods are notsuited for image enhancement in preserving more general structures orother clinically significant structures in reconstructed CT image data.

Despite computationally complex implementation, anisotropic smoothingtechniques have gained attention to reduce noise while preservingstructures of interest. While some prior art anisotropic smoothingtechniques improved in preserving small structures of interest in imagedata, these prior art techniques undesirably require to determine animportance map, which indicates potentially relevant structures thatshould be preserved. Anisotropic diffusion denoising methods and systemsthus remain desired to be improved in preserving clinically significantstructures in reconstructed images using a variety of modalitiesincluding CT and more broadly, multi-dimensional signals.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram illustrating an x-ray computed topographic imagingsystem according to the present invention.

FIG. 2 is a flow chart illustrating exemplary steps involved in theadaptively weighed anisotropic diffusion (AWAD) method or processaccording to the current invention.

FIG. 3A is a diagram illustrating an intensity level corresponding togradient and Laplacian for a certain sharp edge.

FIG. 3B is a diagram illustrating an intensity level corresponding togradient and Laplacian for a certain slow edge.

FIG. 4A is an original reconstructed image before the adaptivelyweighted anisotropic diffusion is applied.

FIG. 4B is a resulted diffused image that has gone through theadaptively weighted anisotropic diffusion (AWAD) according to thecurrent invention.

FIG. 4C is a resulted diffused image that has gone through theadaptively weighted sharp source anisotropic diffusion (AWSSAD)according to the current invention.

DETAILED DESCRIPTION

As described above, the current invention is not limited to CT imagingsystems and is applicable to other imaging systems using modalities suchas MRI, ultrasound and so on. Although the following embodiments of thecurrent invention are described using CT imaging systems, theseembodiments are broadly applicable to multi-dimensional signals.

FIG. 1 is a diagram illustrating an x-ray computed topographic imagingsystem according to the present invention. The projection datameasurement system includes a gantry 1, which accommodates an x-raysource 3 and a two-dimensional array type x-ray detector 5. The x-raysource 3 and the two-dimensional array type x-ray detector 5 arediametrically installed on a rotating ring 2 across a subject, who islaid on a sliding bed 6. The rotating ring 2 rotates around apredetermined rotational axis while the sliding bed 6 moves along therotational axis at a predetermined speed. A gantry/bed controller 9synchronously controls the revolution of the rotating ring 2 and thelinear movement of the sliding bed 6. In the following embodiments,projection data is collected while the x-ray source 3 travels either ahelical trajectory during a helical data acquisition or a circulartrajectory during a circular data acquisition.

The x-ray source 3 emits x-ray flux in a certain configuration such as acone-beam that is directed onto the subject through an x-ray filter 4.An x-ray controller 8 supplies a trigger signal to a high voltagegenerator 7, which applies high voltage to the x-ray source 3 uponreceiving the trigger signal. Furthermore, a system controller 10 exertsan overall control over both the x-ray controller 8 and the gantry/bedcontroller 9 so that x-rays are emitted continuously or intermittentlyat fixed angular intervals from the x-ray source 3 while the rotatingring 2 and the sliding bed 6 are respectively in predetermined rotationand motion.

X-rays that have passed through the subject are detected as electricalsignals by the two-dimensional array type x-ray detector 5. In the x-raydetector 5, a plurality of detector elements is arranged in onedimensional row, and a plurality of the rows is stacked to construct atwo-dimensional array. In certain situation, each detector elementcorresponds to one channel. A data collection unit 11 amplifies theoutput signal from the two-dimensional array type x-ray detector 5 foreach channel and converts to a digital signal so as to produceprojection data.

A processing unit 12 subsequently performs various processing tasks uponthe projection data. For example, the processing unit 12 performs on theprojection data desired operation including but not limited to datasampling and shifting, filtering, backprojection and reconstruction.

The processing unit 12 backprojects the projection data in the imagedomain according to a predetermined algorithm. In general, thebackprojected data is weighted by a distance factor. The distance factoris inversely proportional to the distance L from the x-ray sourceposition to the reconstructed pixel. The distance factor can beproportional to 1/L or 1/L². Also, some additional data redundancyweighting is optionally applied during the backprojection step on thepixel-by-pixel basis. The backprojection step generally obtains datavalue corresponding to the ray through the reconstructed pixel by eitherdata interpolation or data extrapolation. This process can be done in avariety of ways.

The processing unit 12 determines backprojection data reflecting thex-ray absorption in each voxel (three-dimensional pixel) in a 3D imagingreconstruction. In a circular scanning system using a cone-beam ofx-rays, the imaging region (effective field of view) is of cylindricalshape of radius R centered on the axis of revolution. The processingunit 12 defines a plurality of voxels in this imaging region and findsthe backprojection data for each voxel. The three-dimensional image dataor tomographic image data is compiled based upon the backprojectiondata, and the processing unit 12 or other units perform tasks associatedwith denoising of the reconstructed image data. Finally, a displaydevice 14 displays a three-dimensional image or reconstructed tomographyimage according to the three-dimensional image data.

According to the current invention, embodiments of thepost-reconstruction processing unit 13 perform noise reduction in the3-D attenuation map or the reconstructed image that is generated fromthe above described measurements or projection data. In general, theimage quality of reconstructed CT images is expected to improve if noiseis reduced on and around certain structures provided that the resolutionof these structures is at least maintained. Clinically significantstructures in CT images include both large structures and smallstructures. For example, the large structures include bones, the aortalumen and liver while the small structures include small vessels,calcifications and bronchi.

According to the current invention, embodiments of thepost-reconstruction processing unit 13 perform adaptively weightedanisotropic diffusion (AWAD) on the reconstructed images so as to reducenoise while they improve the edge preservation in and around therelevant structures of various sizes. In other words, the noise aroundboth large and small structures is reduced while the edges are betterpreserved according to the embodiment of the current invention.Furthermore, the AWAD has some advantages over the small structures inreducing the noise while preserving the edges. Although many existingprior art anisotropic diffusion techniques can reduce noise in andaround large or relatively slow varying structures, the AWAD bringssubstantial detail of the image back into the diffused image especiallyfor the small structures.

To understand the embodiments as to how to apply adaptively weightedanisotropic diffusion (AWAD) according to the current invention, ageneral concept of anisotropic diffusion is first described with respectto the following equations. Anisotropic diffusion is as an extension ofthe linear diffusion equation by adding a directionally dependentconduction coefficient D as defined in Equation (1a) below:

$\begin{matrix}{\frac{\delta \; u}{\partial t} = {\nabla{\cdot \left( {{D(u)}{\nabla u}} \right)}}} & \left( {1a} \right)\end{matrix}$

where “•” represents inner product of two vectors and D is diffusioncoefficient. A general interpretation of Equation (1a) is that the imageu is diffused over time t in voxels, where the diffusion is steered bythe diffusion tensor D(u) that is based on local image geometry orgradient ∇u. Equation (1a) states that the intensity change per voxel isequal to the divergence of the intensity gradient after a diffusiontensor transformation. As denoted by

$\frac{\delta \; u}{\delta \; t},$

the anisotropic diffusion technique involves that the image u changesover time t in an iterative fashion (i.e. iteration dependent).

In discrete time format, anisotropic diffusion is defined as follows inEquations (1b) and (1c) to explicitly indicate its iterative nature overtime:

$\begin{matrix}{\frac{{u\left( t_{n + 1} \right)} - {u\left( t_{n} \right)}}{t_{n + 1} - t_{n}} = {\nabla{\cdot \left( {D{\nabla{u\left( t_{n} \right)}}} \right)}}} & \left( {1b} \right) \\{{u\left( t_{n + 1} \right)} = {{u\left( t_{n} \right)} + {\left( {t_{n + 1} - t_{n}} \right)\left\lbrack {\nabla{\cdot \left( {D{\nabla{u\left( t_{n} \right)}}} \right)}} \right\rbrack}}} & \left( {1c} \right)\end{matrix}$

where u is the processed or diffused image, and u₀ is the originalimage. u is a function of both position and time. t_(n) and t_(n+1) arerespectively the time instances at iteration n and n+1.

As described above, the anisotropic diffusion technique such as definedin Equation (1) requires that a reconstructed image u is function ofboth space and time. Thus, Equation (2) defines an original image u₀,which in turn defines each voxel at a predetermined initial time t=0.

u(x,y,z,t)_(t=0) =u ₀   (2)

where x, y and z respectively indicates a voxel coordinate in space.Similarly, Equation (3) defines any image u, which in turn defines eachvoxel at a predetermined time t≠0.

u(x,y,z,t)_(t≠0) =u   (3)

where x, y and z respectively indicates a voxel coordinate in space. Inother words, Equation (2) defines a reconstructed image at an initialtime while Equation (3) defines a reconstructed image at a subsequenttime.

Furthermore, Equation (4) defines the Newmann boundary condition whichmeans the derivative in the normal direction of the boundary is 0 andtherefore avoids the energy loss on the image boundary during thediffusion process.

δ_(n)u=0   (4)

The diffusion coefficient D is function of a parameter s as defined inEquations (5a) and (6).

$\begin{matrix}{{D(s)} = \frac{1}{1 + \left( {s/k} \right)^{2}}} & \left( {5a} \right) \\{s = {{\nabla\left( {G_{\sigma}u} \right)}}} & (6)\end{matrix}$

where s is the magnitude of a local gradient of the smoothed u. As notedby G_(σ)u, a reconstructed image u is smoothed by a Gaussian filter Gwith standard deviation σ, k is a threshold or a predetermined value. Ifthe magnitude of the local gradient s at a give pixel is smaller thanthe predetermined threshold value k, the pixel value is treated asnoise. On the other hand, if the magnitude of the local gradient s at agiven pixel is larger than the predetermined threshold value k, thepixel value is treated as a signal for the structures of interest suchas lines or edges.

The predetermined threshold value k is determined based upon thefollowing Equations (5b) and (5c):

k=k₀d   (5b)

Generally, k₀ is set to 2 in Equation (5b) while d is determined asfollows in Equation (5c).

d=average of the gradient magnitude   (5c)

Since the local magnitude of gradient s in Equation (5a) is varying frompixel to pixel, the local magnitude of gradient s is function of spatialpositions. The gradient magnitude defined in (5c) is a single number foran entire reconstructed image.

The above described anisotropic diffusion technique generally workssufficiently well to reduce noise in and around large or relatively slowvarying structures. Unfortunately, the above described anisotropicdiffusion technique alone generally fails to provide detailed edges inand around small or relatively fast varying structures. In other words,the above described anisotropic diffusion technique generally smoothesthe detailed or thin edges in and around the small or relatively fastvarying structures to an extent that the edges are not sufficientlypreserved as clinically significant marks. Consequently, the diffusedimages substantially lose necessary clinical information due to theconventional anisotropic diffusion technique.

To counterbalance the diffusing or smoothing effects on thereconstructed image during the denoising method using an anisotropicdiffusion filter, embodiments of the current invention incorporates anadaptively weighted source term. The weighted source term is adapted foreach instance of iteration and accordingly adjusts the diffused imagefor each instance of iteration. This approach is contrasted to weightingthe finished diffused image at the completion of the iterativeanisotropic diffusion. Consequently, embodiments of a denoising methodand system using an adaptively weighted anisotropic diffusion filteraccording to the current invention bring some of the lost details of theimage back into the diffused image during some instances of iteration ofthe anisotropic diffusion process in order to preserve clinicallysignificant structures in reconstructed images.

The above defined anisotropic diffusion technique is used as ananisotropic diffusion filter or an anisotropic diffusion unit in theembodiments of the current invention. The above defined anisotropicdiffusion is also used in an anisotropic diffusion step in theembodiments of the current invention. In either case, the anisotropicdiffusion is generally implemented in an iterative manner. That is, thediffused image u is evaluated for a current instance of iteration, andthe currently evaluated diffused image u_(t) will be used as a previousdiffused image u_(t−1) in a next instance of iteration.

For each instance of iteration, the embodiments of the anisotropicdiffusion filter or the anisotropic diffusion unit according to thecurrent invention additionally weight the difference between thereconstructed image at the current instance and the originalreconstructed image in an adaptive manner. By the same token, the stepsof the anisotropic diffusion method or process according to the currentinvention also additionally weight the difference between thereconstructed image at the current instance and the originalreconstructed image in an adaptive manner.

The adaptive weight term has been incorporated in the followingEquations (7a), (7b) and (7c) which express the adaptive weightedanisotropic diffusion (AWAD) which is performed in the steps and by theembodiments of the current invention.

$\begin{matrix}{\frac{\partial u}{\partial t} = {{\nabla{\cdot \left( {D{\nabla u}} \right)}} + {W\left( {u_{0} - u} \right)}}} & \left( {7a} \right)\end{matrix}$

where W is an adaptive weight term. In discrete time format, Equation(7a) is optionally defines as in Equations (7b) and (7c) as follows:

$\begin{matrix}{\frac{{u\left( t_{n + 1} \right)} - {u\left( t_{n} \right)}}{t_{n + 1} - t_{n}} = {{\nabla{\cdot \left( {D{\nabla{u\left( t_{n} \right)}}} \right)}} + {W\left( {u_{0} - {u\left( t_{n} \right)}} \right)}}} & \left( {7b} \right) \\{{u\left( t_{n + 1} \right)} = {{u\left( t_{n} \right)} + {\left( {t_{n + 1} - t_{n}} \right)\left\lbrack {{\nabla{\cdot \left( {D{\nabla u}\left( t_{n} \right)} \right)}} + {W\left( {u_{0} - {u\left( t_{n} \right)}} \right)}} \right\rbrack}}} & \left( {7c} \right)\end{matrix}$

The following descriptions of Equations (7a), (7b) and (7c) will befocused on the adaptive weight term since the basic concept ofanisotropic diffusion has been already described with respect toEquations (1a), (1b) and (1c). The adaptive weight term W multiplies thetemporally caused difference of the reconstructed images between theoriginal reconstructed image before any diffusion and the result of thereconstructed image after the current diffusion. In Equation (7a), theabove difference is expressed as (u−u₀). In the discrete time format inEquations (7b) and (7c), the same difference is expressed (u(t_(n))−u₀).

The weight W in Equations (7a), (7b) and (7c) is implemented astime-variant or time-invariant while it is also spatially variant. Inthe following example as defined in Equation (8), the weight W isoptionally illustrated as a time-invariant term. As defined in Equation(8), the weight term W is evaluated as a normalized value between 0 and1.

W=1−e ^(−|L(x,y,z)|)  (8)

where L is a predetermined edge map of the original image u₀. In oneexemplary implementation, Laplacian of Gaussian is used. In otherexemplary implementations, the edge map is produced by other edgedetection methods such as Sobel edge detector (gradient of Gaussian),Curvature edge detector method or a combination of these known edgedetection methods. For example, a Sobel edge detector is given infollowing two filtering masks as defined in Matrixes (8a):

$\begin{matrix}{\begin{matrix}{- 1} & {- 2} & {- 1} \\0 & 0 & 0 \\1 & 2 & 1\end{matrix}\mspace{14mu} {and}\mspace{14mu} \begin{matrix}{- 1} & 0 & 1 \\{- 2} & 0 & 2 \\{- 1} & 0 & 1\end{matrix}} & \left( {8a} \right)\end{matrix}$

In one exemplary implementation using Laplacian of Gaussian as definedin Equations (9) and (10), the edge map based upon Laplacian of smoothedoriginal image Δu, has certain advantage in simplicity and detectingsharper edges and its transit band. Furthermore, the edge map based uponLaplacian of smoothed original image Δu_(s) reduces the slow-changedboundaries and uniform regions to zero-mean or near zero-mean.

L(x, y, z)=Δu _(s)(x, y, z)   (9)

u _(s)(x, y, z)=u ₀(x, y, z){circle around (×)}h(p,q,r)   (10)

where u_(s) is the smoothed version of an original reconstructed imageu₀ by a smoothing kernel h(p, q, r), which is optionally a simplerunning average filter or a Gaussian weighted smoothing filter.

Examples of a simple running average filter or a Gaussian weightedsmoothing filter are formulated in the following Equations (11a) and(11b):

$\begin{matrix}{{{Running}\mspace{14mu} {Average}\text{:}\mspace{14mu} {h\left( {p,q,r} \right)}} = {\frac{1}{z}{\prod{\left( \frac{p}{\tau_{p}} \right){\prod{\left( \frac{q}{\tau_{q}} \right){\prod\left( \frac{p}{\tau_{r}} \right)}}}}}}} & \left( {11a} \right) \\{{{Gaussian}\mspace{14mu} {Smoothing}\text{:}\mspace{14mu} {h\left( {p,q,r} \right)}} = {\frac{1}{Z}{\exp \left( {{- \frac{p^{2}}{2\; \sigma_{p}^{2}}} - \frac{q^{2}}{2\; \sigma_{q}^{2}} - \frac{r^{2}}{2\; \sigma_{r}^{2}}} \right)}}} & \left( {11b} \right)\end{matrix}$

where τ is the length of the running average filter and σ is the sizeparameter for Gaussian filters.

As described above, the anisotropically diffused image is adaptivelyweighted by the weight term which is the product of the normalizedweight value W and the temporally caused difference of the reconstructedimages between the original reconstructed image before any diffusion andthe result of the reconstructed image after the current diffusioninstance in iteration. In other words, the embodiments and the steps ofthe current invention adaptively weight the diffused image at everyinstance of iteration rather than at the completion of the iterativeprocess.

Now referring to FIG. 2, a flow chart illustrates exemplary stepsinvolved in the adaptively weighed anisotropic diffusion (AWAD) methodor process according to the current invention. The following acts ortasks are performed by the above described embodiments including theprocessing unit 12 and the post-reconstruction processing unit 13. TheAWAD method or process does not limit to these steps illustrated in theflow chart.

In one exemplary process of the AWAD, a step S10 generates an edge mapof the original reconstructed image before any iteration of adaptivelyweighed anisotropic diffusion. Some aspects of the edge map generationhave been previously described with respect to Equations (8a), (9),(10), (11a) and (11b). The edge map stores the weight values to beapplied during each instance of iteration of the adaptively weighedanisotropic diffusion. In the above described example, the edge mapgeneration is space-dependent but is time-invariant. That is, althoughthe weight value is determined for a given pixel or a voxel as specifiedby a set of coordinates in the reconstructed image, the weight value isnot varying for the given pixel or voxel over time in one embodiment. Inanother embodiment, the edge map is optionally generated to make theweight both spatially and temporary dependent.

The step S10 also involves some pre-processing of the originalreconstructed image before the edge detection. As described above withrespect to Equations (10), (11a) and (11b), an original reconstructedimage u₀ is smoothed by a smoothing kernel such as h(p, q, r).

Lastly, the step S10 generates an edge map of certain characteristics.While the generated edge map provides the weight values in a weightingprocess, certain characteristics of a particular edge detector generallydetermine the weight values to be applied in the weighting step. Forthis reason, a selected edge detector partially affects an ultimateoutcome of the adaptively weighed anisotropic diffusion processaccording to the current invention. The edge detectors include a varietyof methods such as gradient map, Laplacian, curvature, Eigen values ofstructural tensors and Hessian matrices. The edge detector isalternatively implemented based upon a combination of the abovedetectors.

Still referring to FIG. 2, a step S20 performs a predeterminedanisotropic diffusion process on the reconstructed image data such as aseed image according to the current invention. An exemplary step of theanisotropic diffusion process has been previously described with respectto Equations (1a), (1b), (1c), (2), (3), (4), (5a), (5b), (5c) and (6).

In a step S30, the anisotropic diffusion process weights the diffusedresult in the step S20 according to the current invention. The adaptiveweighting is accomplished by tasks associated with Equations (7a), (7b),(7c) and (8). Equations (7a), (7b), (7c) includes terms associated withthe anisotropic diffusion process itself and the adaptive weightingterms. Although the steps S20 and S30 are illustrated as two sequentialsteps in the exemplary flow chart, the adaptively weighted anisotropicdiffusion process is implemented in a variety of ways according thecurrent invention. In an alternative embodiment or process, the stepsS20 and S30 are optionally implemented as simultaneous steps that may beperformed in parallel. In any case, regardless of a particularimplementation, the diffused result in the step S20 is weightedaccording to the adaptive weighting in the step S30 to generate anadaptively weighted result for each instance of the iteration. Equation(8) defines the normalized weight value W based upon the edge map thathas been generated in the step S10.

As mentioned before with respect to the step S30, the adaptive weightedanisotropic diffusion according to the current invention is iterativelyperformed in one preferred process. In this regard, a step S40determines as to whether or not the adaptive weighted anisotropicdiffusion in the steps S20 and S30 should be terminated. If the step S40determines that the adaptive weighted anisotropic diffusion in the stepsS20 and S30 should be terminated, the illustrative steps end itsprocess. On the other hand, if the step S40 determines that the adaptiveweighted anisotropic diffusion in the steps S20 and S30 should not beterminated, the illustrative steps continue its process by iterativelyperforming the steps S20 and S30. As described with respect to Equations(7a), (7b) and (7c), the reconstructed image at a subsequent timeu(t_(n+1)) is now assigned to the reconstructed image of a previous timeu(t_(n)) in the current instance of the iteration before a next instanceof the iteration is performed.

In an alternative implementation, the step S30 has an optional step ofapplying a predetermined sharp source filter according to the currentinvention. Equations (12) and (13) define a process of adaptive weightedsharp source anisotropic diffusion (AWSSAD) according to the currentinvention.

$\begin{matrix}{\frac{\partial u}{\partial t} = {{\nabla{\cdot \left( {D{\nabla u}} \right)}} + {W\left( {{S\left( u_{0} \right)} - u} \right)}}} & (12)\end{matrix}$

where S(u₀) is the sharpened version u₀ using a predetermined unsharpfilter as defined in Equation (13).

$\begin{matrix}{{S\left( u_{0} \right)} = \left( {{\frac{C}{{2\; C} - 1}u_{0}} - {\frac{1 - C}{{2\; C} - 1}u}} \right)} & (13)\end{matrix}$

w is a parameter of unsharp filter. At C=0.8, reasonably good resultshave been obtained.

The step S40 determines as to whether or not the adaptively weightedanisotropic diffusion should be further iteratively performed based upona predetermined set of termination conditions. Commonly, there are twoways to determine the number of iterations or the conditions fortermination: One termination is based upon comparing themean-square-errors between current iteration result and previousiteration result. If the difference is smaller than a preset terminationvalue, the processing stops. The other termination is based upon apredetermined noise assessment tool to determine a noise level at eachinstance of the iteration. If the current noise level is reduced by apreset noise reduction goal such as a 50% noise level from the original,the iteration stops.

In one implementation of the adaptively weighted anisotropic diffusionprocess, if iterative methods are used to calculate an image at aspecific scale t in the scale-space by choosing a specific Δt and byperforming t/Δt iterations to calculate the diffused image resultu(t_(n+1)) based upon the previous diffused image result u(t_(n)) from aprevious iteration, the step size Δt is optionally determined accordingthe following Equation (14) for stability.

Δt<1/(2N)   (14)

where N is the number of dimensions of the reconstructed image. Basedupon Equation (14), the value of Δt, which limits the upper bound of thestep size. If the step size Δt is too small, the adaptive weightedanisotropic diffusion process substantially slows down the overalldenoising process. In one implementation, the step size Δt is setapproximately between ¼˜1 of the upper bound. For example, the step sizeΔt takes a value such as ⅛, 1/12 or 1/16 for N=3 in denoising the 3Dreconstructed images.

Now referring to FIGS. 3A and 3B, diagrams illustrate an intensity levelcorresponding to gradient and Laplacian for sharp and slow edges. Asdescribed above with respect to the step S10, an edge map is generatedwith certain characteristics according to a particular edge detector. Inthis regard, a selected edge detector partially affects an ultimateoutcome of the adaptively weighed anisotropic diffusion processaccording to the current invention. In one implementation, an edge mapis generated to emphasize small or sharp edges rather than thick orslow-varying edges or boundaries. One advantage of the adaptivelyweighted anisotropic diffusion is to retain small edges while the imagesare denoised. For this purpose, Laplacian is used in generating an edgemap.

FIG. 3A, a diagram illustrates an intensity level corresponding togradient and Laplacian for a certain sharp edge. A top graph illustratesan intensity level of an exemplary sharp edge. The intensity changesquickly across the zero intensity level. A middle graph illustrates agradient level that corresponds to the top exemplary sharp edge. Thegradient is at the highest level when the edge changes its intensitylevel from a negative value to a positive value at zero. A bottom graphillustrates an intensity level of Laplacian, which is a second orderderivative in one dimensional case. For a rising edge or a falling edgeof the exemplary sharp edge, a first order derivative gives either apositive upward bump or a negative downward bump. The peak or valley ofthe bump is the place where the first order derivative is at maximum inabsolute value. The peak or valley of the first order derivative isactually the zeroes in second order derivative, and at two sides of thezeros, there are two bumps (one positive and one negative) in secondorder derivative or Laplacian. In one implementation, an absolute valueof the second derivative is taken as the weight at edges so that moreweight is given to the pixels within the band along the edges topreserve the edge information. Although the Laplacian weight at edge iszero, the edge pixels are substantially preserved by the anisotropicdiffusion process. Only sharp edges have larger Laplacian values, andthese edges are substantially preserved by using Laplacian weightvalues.

FIG. 3B, a diagram illustrates an intensity level corresponding togradient and Laplacian for a certain slow edge. A top graph illustratesan intensity level of an exemplary slow edge. The intensity changesslowly across the zero intensity level in comparison to the top graph inFIG. 3A. A middle graph illustrates a gradient level that corresponds tothe top exemplary slow edge. The gradient also changes slowly and with aless intensity level when the edge changes its intensity level from anegative value to a positive value at zero. A bottom graph illustratesan intensity level of Laplacian, which shows almost no change in theintensity level. In contrast to a sharp edge, slow varying edges orlarge edges having thick boundaries generally have small Laplacianvalues compared to the noise. For this reason, Laplacian issubstantially ignored in one implementation.

FIGS. 4A, 4B and 4C illustrate some effects of the adaptively weightedanisotropic diffusion on an original image according to the currentinvention. FIG. 4A is an original reconstructed image before theadaptively weighted anisotropic diffusion is applied. FIG. 4B shows aresulted diffused image that has gone through the adaptively weightedanisotropic diffusion (AWAD) according to the current invention using athreshold coefficient=2 with number of iterations=4 and step size=⅛.FIG. 4C shows a resulted diffused image that has gone through theadaptively weighted sharp source anisotropic diffusion (AWSSAD)according to the current invention using threshold coefficient=2 withnumber of iterations=13, step size=⅛ and w=0.8.

In summary, the embodiments of the current invention apply a negativefeedback mechanism to adjust the edge and texture features in areconstructed image based upon the information in the original image sothat the final result is a balanced, converged and featured image at areduced level of noise.

While certain embodiments have been described, these embodiments havebeen presented by way of example only, and are not intended to limit thescope of the inventions. Indeed, the novel methods and systems describedherein may be embodied in a variety of other forms; furthermore, variousomissions, substitutions and changes in the form of the methods andsystems described herein may be made without departing from the spiritof the inventions. The accompanying claims and their equivalents areintended to cover such forms or modifications as would fall within thescope of the inventions.

What is claimed is:
 1. A medical imaging apparatus comprising: ananisotropic diffusion unit configured to execute anisotropic diffusionfor medical image data comprising a plurality of signal values, eachsignal value corresponding to coordinates in a space; and a weightgeneration unit configured to generate a plurality of weights, whereinthe weight for each signal value depends on an edge value correspondingto the coordinates of the respective signal value, wherein theanisotropic diffusion unit repeatedly executes the anisotropic diffusionfor the medical image data by using the weights by a predeterminednumber of iterations, wherein the anisotropic diffusion unit is furtherconfigured to execute the anisotropic diffusion using inhomogeneousanisotropic nonlinear filter, and an inhomogeneous term in theinhomogeneous anisotropic nonlinear filter is a term that multiplies adifference between reference medical image data and medical image datasmoothed by the anisotropic diffusion by the weights respectivelycorresponding to the signal values.
 2. The apparatus of claim 1, whereinthe inhomogeneous anisotropic nonlinear filter is given byu(t _(n+1))=u(t _(n))+(t _(n+1) −t _(n))[∇•(D∇u(t _(n)))+W(u(t _(n))−u₀)] where u(t_(n+1)) is a signal value at the (n+1)th iteration,u(t_(n)) is a signal value at the nth iteration, u₀ is an initial signalvalue before application of the inhomogeneous anisotropic nonlinearfilter, and W is a weight corresponding to the initial signal value. 3.The apparatus of claim 1, wherein the medical image data comprisesvolume data, and the signal value comprises a voxel value.
 4. Theapparatus of claim 1, wherein the reference medical image data comprisesmedical image data before application of the anisotropic diffusion. 5.The apparatus of claim 1, wherein the weight generation unit is furtherconfigured to generate the weights by using Laplacian edge detection. 6.The apparatus of claim 1, wherein the weight generation unit is furtherconfigured to generate the weights by using Sobel edge detection.
 7. Theapparatus of claim 1, wherein the weight generation unit is furtherconfigured to generate the weights by using curvature edge detection. 8.The apparatus of claim 1, wherein the weight is invariant in theanisotropic diffusion repeated by the predetermined number ofiterations.
 9. The apparatus of claim 1, wherein the medical image datacomprises medical image data smoothed in advance before the anisotropicdiffusion.
 10. The apparatus of claim 1, wherein the reference medicalimage data comprises medical image data sharpened based upon medicalimage data before the anisotropic diffusion and a predeterminedparameter.
 11. The apparatus of claim 10, wherein the sharpened medicalimage data is given by${S\left( u_{0} \right)} = \left( {{\frac{w}{{2\; w} - 1}u_{0}} - {\frac{1 - w}{{2\; w} - 1}u}} \right)$where S(u₀) is a signal value of the sharpened medical image data, u₀ isan initial signal value before application of the inhomogeneousanisotropic nonlinear filter, u is a signal value at the nth iteration,and w is the predetermined parameter.
 12. The apparatus of claim 1,wherein the anisotropic diffusion unit is further configured to store apredetermined time scale, the anisotropic diffusion is executed at aspecific time interval, and the predetermined number is a numberobtained by multiplying the time scale by the specific time interval.13. The apparatus of claim 12, wherein the specific time interval issmaller than a reciprocal of twice the number of dimensions of themedical image data.
 14. A medical image processing method comprising:executing anisotropic diffusion for medical image data comprising aplurality of signal values, each signal value corresponding tocoordinates in a space; and generating a plurality of weights, whereinthe weight for each signal value depends on an edge value correspondingto the coordinates of the respective signal value, wherein theanisotropic diffusion is repeatedly executed for the medical image databy using the weights by a predetermined number of iterations, usinginhomogeneous anisotropic nonlinear filter, and an inhomogeneous term inthe inhomogeneous anisotropic nonlinear filter is a term that multipliesa difference between reference medical image data and medical image datasmoothed by the anisotropic diffusion by the weights respectivelycorresponding to the signal values.